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ABSTRACT 

Many objects studied in astronomy follow a power law distribution function, for example 
the masses of stars or star clusters. A still used method by which such data is analysed is to 
generate a histogram and fit a straight line to it. The parameters obtained in this way can be 
severely biased, and the properties of the underlying distribution function, such as its shape or 
a possible upper limit, are difficult to extract. In this work we review techniques available in 
the literature and present newly developed (effectively) bias-free estimators for the exponent 
and the upper limit. Furthermore we discuss various graphical representations of the data and 
powerful goodness-of-fit tests to assess the validity of a power law for describing the distribu- 
tion of data. As an example, we apply the presented methods to the data set of massive stars 
in R136 and the young star clusters in the Large Magellanic Cloud. For R136 we confirm the 



result of Koen| ( |2006| l of a truncated power law with a bias-free estimate for the exponent of 



2.20 ± 0.78 / 2.87 ± 0.98 (where the Salpeter-Massey value is 2.35) and for the upper limit of 
143 ± 9M / 163 ± 9M , depending on the stellar models used. The star clusters in the Large 
Magellanic Cloud (with ages up to 10 7 5 yr) follow a truncated power law distribution with 
exponent 1 .62 ± 0.06 and upper limit 68 ± 12 x 1O 3 M . Using the graphical data representa- 
tion, a significant change in the form of the mass function below 10 2 5 M can be detected, 
which is likely caused by incompleteness in the data. 

Key words: methods: statistical - methods: data analysis - stars: luminosity function, mass 
function - galaxies: star clusters 



1 INTRODUCTION 

Many astronomical objects are distributed according to a power 
law. The probably most prominent example is the mass function 
of stars more massive than 0.5 M a with the Salpeter-Massey expo- 
nent of 2.35. Further examples are the mass functions of young star 
clusters and of molecular clouds. Modern observational techniques 
and state-of-the-art models provide data such as stellar masses with 
unprecedented accuracy. However, the statistical analysis of those 
data is not yet always optimal. The technique of binning the data 
suffers from losing a lot of information. The grouping of data into 
cells instead of using every data point obscures details of the ob- 
served distribution. This is an especially serious problem in the 
upper range, where the bins are only sparsely filled. Furthermore 
the obtained estimates of the slo pe can be severely biased (see e.g. 
Maiz Apellaniz & Ubeda 20051. A method based on a particular 
graphical display of the data which avoids grouping and allows one 
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an estimate of the upper limit was given by lKoenlprjOo) . Another 
successful approach is to use the Maximum Likelihood method, 
which has been applied by | Jauncey| ( |1967| > on extragalactic radio 
sources. ICrawford et"aTT1 < | 1970| > derived a Maximum Likelihood es- 
timator for the exponent without grouping the data and including 
an upper limit. 

A further step in data analysis, equally important as estimat- 
ing the parameters, is the validation of the assumed power law form 
of the distribution. The simplest way to do this is to look at the 
histogram of the data in a double logarithmic plot. If this plot ap- 
pears to be linear then the consistency of the data with a power 
law is concluded. But the significance of deviations from linearity 
are hard to state in an objective way by mere visual inspection. A 
further, more elaborate way is to apply a goodness-of-fit test such 
as the Kolmogorov-Smirnov test. If the calculated test statistic lies 
in some acceptance range then also consistency is concluded. But 
it is possible that the test statistic calculated with data stemming 
from an alternative hypothesised distribution similar to the power 
law might as well fall in the acceptance range. The test then fails 
to produce the right result since it has not enough "power" to dis- 
criminate. Therefore the "power" properties of a goodness-of-fit 
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test have likewise to be examined. Such a study is - to our knowl- 
edge - not yet available in the astronomical literature. 

In this work we describe estimation methods and compare 
their biases and variances (Section |3J. Since the data may stem 
from a truncated power law we focus on estimators which can be 
used in this case. In the second part (Section |4| we investigate the 
question whether the data are consistent with the assumed power 
law distribution. As informal aids to answer this we discuss var- 
ious plotting recipes (Section |4.1[ ). For an objective decision we 
present goodness-of-fit tests with a study of their discrimination 
power (statistical power) under the hypotheses of a truncated and 
infinite power law. Finally, in Section [5] we will apply the intro- 
duced methods on the massive stars in R136 and the young star 
clusters in the Large Magellanic Cloud. 




log 10 x zmax 

Figure 1. Complementary cumulative DF (CCDF) plot for an infinite (dot- 
ted line) and truncated (solid line) power-law pdf (a = 2.35, *min = 1. 
^MAX = 150, shown by the vertical thick bar). For the truncated case a char- 
acteristic turn-down appears at the upper end. 



2 GENERAL RESULTS, DEFINITIONS AND NOTATION 
2.1 The power law distribution 

In this work we only consider power law distribution functions 
(DF) with a negative exponent, — a (a > 1). By convention the sign 
is separated from the absolute value. Besides the exponent such a 
distribution is further parametrised by the lower and upper limit, 
*min an d *max- The probability density is then given by 



p(x;a,x Mm ,x MAX ) 
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and the cumulative distribution function (DF) is 



P(x) 
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The family of distributions given by eq. ^ includes the "in- 
finite" or "not truncated" distributions with infinite upper limit, 
p(x;xmax = °°) P~(jc), which is also known in the (non- 
astronomical) literature as the Pareto distribution. The density func- 
tion reads then 



X MIN 

and the cumulative distribution is 

l-a 
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An useful property of the power law distribution is its rela- 
tion to the exponential distribution. By a logarithmic transforma- 
tion the power law distribution becomes proportional to an expo- 
nential, x~ a = e~ alo S«*. Due to this proportionality it is possible 
that some techniques for estimation and testing, which were devel- 
oped for the exponential distribution, can be used for the power law 
distribution. 



3 ESTIMATING THE PARAMETERS 

There exist in the literature a variety of methods to estimate the 
parameters, exponent and upper limit, of a power law distribution. 
In this first part of the paper we describe them and compare their 
properties. 



3.1 Methods 

3.1.1 Binning 

A particularly simple method is to fit a linear relation to the data 
gr ouped in bins of constant size in logarithmic space. As shown 
by Mafz Apellaniz & Ubeda ( 2005 1 this method can yield biased 
results, i.e. results which systematically deviate from the actual pa- 
rameter, and do not allow one to estimate a possible upper limit. 
A solution to avoid biased results was given by Mafz Apellaniz 



&Ubedal|2005b, which modified the binning scheme to a constant 



number of data points per bin and fitted the expected number in- 
stead of performing a linear regression. In this way the estimate for 
the exponent, S, can be obtained, together with an estimate of its 
uncertainty, which is consistent with the sampling variance of S. 
In extension of their work we investigate the properties of the es- 
timate for the upper limit, derived from the normalisation constant 
of the frequency distribution, k = n( 1 — a) / (*max — x min ) > which 
is given by 



*MAX 
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-X, 
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(5) 



where the smallest observation (Xn\) is used as an estimate for 

*MIN- 

Not only the choice of constant-size or variable-size bins has 
influence on the results of binning, but also the number of bins. 
D' Agostino & Stephens ( 1986 1 give as the optimal number of bins 
2n 1 ^ for n data points. A smaller number of bins reduces the bias 
but increases the standar d deviation (cf. Table 1 to 3 of |Mafz| 
Apellaniz & Ubeda|2005 i. 



3.1.2 Complementary cumulative distribution function plot 



|Koen| ( (2006[ > presented a method to estimate both the exponent and 
the limits of a power law. This method is based on a particular 
graphical representation of the data, the complementary cumulative 
DF (CCDF) plot (Fig.[TJ. Data stemming from an infinite power law 
follow a linear relation with slope 1 — a in a plot of log ( 1 — P*, (x)) 
versus logx, as can be seen easily by taking the log of eq.|4] For a 
truncated power law, a turn-down appears at the high end. Estimates 
for the exponent and limits are obtained by fitting log(l — P(X/ t \)) 
(with P(x) from eq.^and the ordered data X^) to log(l - *=jp), 
with using '~^ ,s for the empirical cumulative distribution function. 
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3.1.3 Beg 's estimator 

A power-law distribution is closely related to the exponential dis- 
tribution. Therefore it is possible to apply the uniformly minimum 
variance unbiased estimators for the slope and limits of a truncated 
exponential distribution, developed by Beg] (|1982| [T9 83 ) to log- 
transformed power law data, as shown by Beg (1983). Although 
these estimators are theoretically an optimal solution, they are only 
partially practicable, since their computation is numerically diffi- 
cult and impossible for large data sets. We therefore developed a 
recursive form which is applicable to arbitrarily large data sets. The 
original and recursive formulae are given in the Appendix. 




3.1.4 Maximum Likelihood (ML) estimator 

The Maximum Likelihood estimator for the exponent was given by 
|Crawford et al.| ( 19701, who also included an upper limit. In our 
case the upper limit is intrinsic to the distribution function, i.e. the 
upper limit needs to be estimated simultaneously with the expo- 
nent. The likelihood function for a random sample of size n from a 
truncated power law DF is 



Qpfx;; a;x M IN,^MAx) 
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The estimator for the exponent is obtained by maximising the log- 
likelihood, 



logJzf = filog(l — a)— iilog(.T 
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The maximisation can be performed by finding the root of the 
derivative with respect to a of eq. [8] The estimator CLml is then 
the solution of 
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with Y = minX;, Z = maxXj and T = logX,. 

The ML estimates for the upper limit Xmax = vaaxX, (see 
e.g. |Aban et al.|200"6| >. It is obvious that this estimate will be biased 
since the upper limit is larger than the largest data point. 

3.1.5 Bias-free estimators based on the maximum likelihood 
estimator 

It is possible to construct a minimum variance unbiased estimate of 
the exponent from the maximum likelihood estimate, as shown by 
|Crawford et al.H1970) or |Baxter|(l98(ty . For the infinite case the 
ML estimator for the exponent is given by 



a — 1 



(10) 



r-nlog e F 

with Y = minX, or the given lower limit, and T = £™ =1 log^X,-. The 
unbiased estimator is then 

S'-l = — (S-l) (11) 
n 

(if both the exponent and the lower limit should be estimated then 



(n - 2) jn has to be used i Baxter 1980 1). 

The simple relation between the ML estimator and the unbi- 
ased estimator for an infinite power law suggests a similar relation 
for the truncated case. However, for a truncated power law a closed 



Figure 2. Distributions of the modified ML estimates for the exponent 
(histogram), derived from Monte-Carlo samples of size 1000 for three input 
values (1.5, 2.0, and 2.5). They follow a Gaussian with mean and variance 
derived from the samples. For larger exponents the variance increases. 



form of the ML estimator is not available. This makes a proof of an 
unbiased estimator very difficult and maybe even impossible, since 
the distribution of the estimate cannot be calculated analytically. 
Nevertheless, it is not unreasonable to assume that a simple modifi- 
cation of the ML estimate also leads to unbiased results. A different 
pre-factor, depending only on the number of data, should give the 
expected result. We found that 



-(%£-!) 



(12) 



(MML = Modified Maximum Lik elihood) provid es quasi-bias-free 

(n — 3 be- 



Baxter 



1980) or 



n-3 



estimates. The pre-factors ^= of 
cause there is an additional parameter, the upper limit) lead to bi- 
ased results. The distribution of the exponents estimated using this 
method follows a Gaussian, as can be seen in Fig. [2] with an in- 
creasing variance for an increasing exponent. 

The bias of the ML estimate for the upper limit can also be sig- 
nificantly reduced by appropriate modifications. Hannon & Dahiya 
( 1999 1 developed such a modified estimator for the exponential dis- 
tribution. This estimator can also be used for the power law distri- 
bution and takes then the form (with the ML estimate of the expo- 
nent replaced with the bias-reduced form) 



*MAX 



with 



1 + - 



(l-«MML)log ? 



(13) 



(14) 



where Xn\ is the smallest and Xi n \ 



is the largest data point. The 
properties of the modified estimate are discussed in the next Sec- 
tion. 



3.2 Performance of the estimators 

After introducing a number of methods of estimation, we compare 
their properties. The quality and usability of an estimator is deter- 
mined by several factors. A main demand is that an estimate is on 
average equal to the actual parameter, i.e. bias-free. Also, the vari- 
ance of the estimate should be as small as possible and it should be 
numerically robust. 

To study these properties we carried out a set of Monte-Carlo 
experiments, each of size 1000, with parameters in the typical range 
of astronomical applications. The values for the exponent range 
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Number of data Number of data 



Figure 3. In the left panel the results for the exponent are shown, on top the average bias (calculated using eq. 1 1 5| the horizontal lines mark ±0.025) and 
below the average standard deviation (eq. |16} of the estimated exponents. In the right panel, average estimates of the upper limit (lines mark the true values) 
are displayed in the upper part and below the average relative bias is shown (eq. |17[ . The parameter combinations are given in the text, Sec. |3.2| The symbols 
refer to: A constant size binning using linear regression; T variable size binning using x 2 ', ♦ CCDF plot fitting. ML estimator without including truncation; 
• Beg's estimator; (This estimator starts to fail for n ^ 150, dots below the x-axis indicate a failed experiment); o Beg's estimator in the recursive form; x ML 
estimator; i( Modified ML estimator; 



from 1.6 to 2.85 in steps of 0.25. For each exponent four pairs 
of limits were used ({0.5, 150} and {10, 150} corresponding to the 
stellar mass function, { 10 3 , 10 5 } and {10 4 , 10 6 } corresponding to 
the mass function of young star clusters). The last varied parameter 
was the number of data (50, 100, 300). For the binning methods the 
numbe r of bins was chosen according to D'Agostino & Stephens 
1 1986 i (2n 2 / 5 ), which gave 9, 12 and 19 bins, respectively. 

As diagnostics for the performance in estimating the exponent 
we choose the average bias of an estimator for a given parameter 
set, 



B(a) 



1 



1000 



1000 



E (%-«)> 



i=1 



and the standard deviation, 
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1000 
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(15) 



(16) 



The left panel of Figure [3] summarises the results for the bias. 
Two horizontal lines at ±0.025 embrace the region, in which we 
consider the bias as negligible. The general trend is that the bias 
decreases with an increasing number of data (except for the ML 
estimator wich does not include a truncation). The corresponding 
results for the standard deviation also decrease with larger size of 
the data set. 

Variable-size binning gives effectively bias-free exponents for 
samples having a moderate size or larger. The method of Koen is bi- 
ased towards lower exponents. The results from Beg's estimator are 
very good, but the method fails for large data sets, even in the recur- 
sive form. A maximum likelihood estimate without considering the 
truncation can lead to a significantly overestimated exponent. But 
when the truncation is considered, the bias is small and effectively 
vanishes if our modified version is used. 



The results for the estimates of the upper limit are shown in the 
right part of Fig. [3] Generally it can be observed that a larger upper 
limit leads to larger absolute deviations in the estimate. Because 
the upper limits used in this study span a wide range of values it is 
not convenient to compare the absolute biases as for the exponents. 
The relative bias (also displayed in Fig.|3|, 

1000 

E 



S( x MAx) 
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1000 



*MAX.i — *MAX 
*MAX 



(17) 



is a better measure of trends. Furthermore the normalised distribu- 
tions of the estimates are shown in Fig. [4] for two parameter sets 
(a = 2.0, n = 100, limits {10, 150} and {1000, 100000}). The his- 
togram and a Gaussian with mean and variance (cr) calculated from 
the Monte-Carlo sample are rescaled by x' = Jt ~^ 1AX and the y-axis 
is scaled such that the peak of the Gaussian is 1 . If the estimator is 
not biased and can be approximated by a Gaussian, then the nor- 
malised distribution should follow a Gaussian with zero mean and 
unit variance. 

The upper limit is underestimated by using the normalisation 
constant of the variable-size binning method. The results of fitting 
the CCDF plot are peaked around the input value, but can have a 
long tail of very high estimates (for the limits {1000, 100000}). If 
in the CCDF plot the data show no strong curvature at the upper 
end, then the estimated upper limit is very large. The distribution 
of estimates obtained with Beg's estimator are in reasonable agree- 
ment with a Gaussian, but not completely symmetric around the 
mean. If the largest data point is used (i.e. the direct ML estimate), 
then the upper limit is underestimated, with a distribution limited 
by the actual value. With the modification (eq. 1 13^ the distribution 
becomes similar to the one of Beg's estimator, spreading around the 
true value. Although not completely symmetric around the mean it 
can be sufficiently approximated by a Gaussian, and has no outliers 
as fitting the CCDF plot. 
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Figure 4. Distribution of the estimated upper limits for the different meth- 
ods (histograms; parameters: a = 2.0; n = 100; solid: limits {10,150}, 
dashed: limits { 1000, 100000}.) The x-axis has been scaled such that a un- 
biased estimate should follow a Gaussian of zero mean and unit variance 
(dotted). Also shown is a Gaussian with mean and variance derived from 
the estimates. 



In summary, when both the exponent and the upper limit of a 
truncated power law should be estimated, our modified ML method 
performs best in terms of bias and stability, being similar to Beg's 
uniformly minimum variance unbiased estimator, but without the 
numerical instability, gives the best results in terms of bias and sta- 
bility. 



4 IS A POWER LAW CONSISTENT WITH THE DATA? 

For a thorough analysis of data which are assumed to stem from a 
power-law distribution it is not sufficient just to estimate the pa- 
rameters. The parameter estimation answers the question which 
power-law fits the data best, but leaves open whether the data are 
originating from a power law at all. Or to put it differently: is the 
(truncated) power law a good parent distribution function of the 



data? The need to answer this question has already been stressed by 
|Crawford et al.| l|l970 1. This question can be addressed by a graph- 
ical inspection of the data, which is discussed in the next Section. 
After the informal visual methods more objective goodness-of-fit 
techniques are discussed. 



4.1 Graphical inspection of the data 

A common approach to find the parent DF of a data set is to use 
a histogram as a non-parametric estimate of the form of the parent 
DF. If in a logarithmic plot the histogram of e.g. stellar masses is 
a straight line, a power-law is usually assumed as the parent DF. 
However, a power-law is a heavy-tailed distribution and has only 
a few counts per bin in the tail. Thus the scatter in a histogram 
is large in the upper regime and makes deviations from a power- 
law hard to detect. Alternative, heavy-tailed distributions lead to 
nearly indistinguishable histograms. It is for example not possible 
to decide whether the power-law is truncated or not. Therefore a 
histogram only allows us to roughly determine the parent DF. 

A display of the data which avoids grouping them into cells is 
the probability-probability (or percentile-percenile, PP) plot. For 
the PP plot the data have first to be sorted in ascending order, 
Xf(\ < Xr i+ iy The x-values then follow as the "theoretical" per- 
centiles, which are the values of the cumulative DF for the i-th 
data point, P(X^,oc,xmin,xmax), calculated with the estimated 
parameters. As y- values the "empirical" percentiles are used, given 



J by 



1-0.5 



Both axes range from to 1, and when the data lie on 



the diagonal then they agree with the null-hypothesis. Hypothe- 
ses other than the null hypothesis can be shown by plotting the 
pairs {x=al tentative cumulative DF,y=null cumulative DF}. An ex- 
ample for a PP plot with an infinite power law as null hypothesis 
(diagonal) is shown in Fig. [5] but with data generated from a trun- 
cated power law. The curve for the alternative hypothesis of a trun- 
cated power law (of the same exponent, solid line) barely deviates 
from the diagonal and does not allows one to distinguish the in- 
finite and truncated versions. In this plot the acceptance region of 
the Kolmogorov-Smirnov (KS) test can be directly shown as paral- 
lels to the diagonal, in Fig.|5]calculated with a significance level of 
5%. The data do not exceed this region, not even in the tails, giving 
evidence for the known insensitivity in the tails of the Kolmogorov- 
Smirnov statistic. 

Before showing a way to improve the insensitivity of the PP 
plot in the tail we shortly compare it with other possible plots which 
show all available data (see e.g. |Chambers et aT] ( |1983| > or |Wilk &| 
Gnanadesikan j!968| ). A plot which uses e.g. the X^ as x-values 
(the "empirical cumulative density plot") has a curved reference 
line for both an infinite and truncated power law, with only small 
and not well perceptible differences between them. Following the 
suggestion of Koen (2006), in a plot of log(l —P(X)) against log A" 



(the CCDF plot, see Fig. [T] and Sec. |3.1.2| l the data should only 
show a curvature for a truncated power law. But from the results 
for the estimator based on this plot, the scatter in the highest data 
points can be large, and a graphical goodness-of-fit criterion would 
not be very sensitive. A third alternative to the PP plot would be 
a plot of the inverse cumulative DF (P~ 1 ( '~^' 5 ), giving the ex- 
pected value for the data point X^\) against the ordered data X/ t \ 
("quantile-quantile" plot). This plot has a linear reference, the up- 
per tail would be curved when constructed for the infinite null hy- 
pothesis but with truncated data. Again, as for the complementary 
cumulative density plot there is large scatter, and additionally the 
inverse cumulative DF has to be calculated. 
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P (X (a ) or P(x) S (^(%)) or s 

Figure 5. Example for a percentile-percentile plot (PP, left) and a stabilised percentile-percentile plot (SPP, right) for the null hypothesis of an infinite power 
law (=diagonal), using 100 data points sampled from a truncated power law (a = 2.35, .vmin = 1 and jcmax = 150). Also shown are the curves for a truncated 
power law (solid line, parameters as estimated, 2 = 2.35 and £max = 149). The acceptance region of the (in the right plot stabilised) Kolmogorov-Smirnov 
statistic (significance level 5%) is given by the two parallels to the diagonal. The data lie within this region in the PP plot, wherefore from the PP plot the 
infinite power law cannot be significantly rejected. After stabilisation the KS test is more powerful and thus allows us to detect truncation in contrast to the PP 
plot. 



4.2 The stabilising transformation and the SPP plot 

Goodness-of-fit methods based on the empirical cumulative DF, 
such as the PP plot or the KS test, have the advantage that their in- 
trinsic properties do not depend on the actual choice of the null hy- 
pothesis. A PP-plot for e.g. Gaussian variates looks identical (mod- 
ulo the random scatter) to a PP plot for power-law data. The reason 
for this is that by taking the cumulative DF the data are transformed 
to uniformly distributed variates, if they are following the null hy- 
pothesis. Therefore the location and variances of the points in the 
PP plot are independent of the null hypothesis distribution, as is 
the distribution of the KS statistic. This transformation reduces the 
goodness-of-fit task from arbitrary DFs to testing for uniformity. 
However, the variances of uniform ordered variates are not inde- 
pendent from their position: the scatter in the PP plot is larger in the 
middle than in the tails. Thus a test which measures the differences 
between the expectation and the empirical values of uniformly dis- 
tributed ordered data will be dominated by the points with the larger 
variances. Hence the insensitivity of the KS-test in the tails. 

A way to overcome the unequal variances was introduced by 
Michael ( 1983 1. The stabilising transformation of uniform variates 
u (= the cumulative DF), 

2 

Sq(u) = — arcsin(v / «) (18) 

gives asymptotically equal variances of the transformed ordered 
variates. In a stabilised PP (SPP) plot every part of the plot has the 
same weight and no region is particularly emphasised. Although 
the distribution of the Sq{u) is not uniform any more, tests based 
on the differences between expectation and empirical value can still 
be used (as long as they do not use other properties of the uniform 
distribution). These transformed tests are equally sensitive to every 
part of the distribution function. 

However, for testing the tail-behaviour of a DF it is useful to 



emphasise the tail. This can be achieved by using only a half trans- 
formation, which is possible because So is symmetric around the 
point {0.5,0.5} and the interval [0.5, 1] is mapped onto [0.5, 1]. A 
one-sided transformation of the percentiles to stabilise a right-tailed 
distribution consists then of three steps. First, map the interval [0, 1] 
on [0.5, 1], then use So, and lastly map [0.5, 1] back on [0, 1]. The 
formula for this is 

S(u) = 2S (0.5 + 0.5m)- 1. (19) 

We use S instead of So in the SPP plot (Fig. [5} and in the goodness- 
of-fit tests which are related to it, because of the one-tailed power 
law distribution. 

4.3 Goodness-of-fit Tests 

A goodness-of-fit test provides an objective way to "measure" the 
agreement of the fit with the data. We follow here the Neyman- 
Pearson ansatz of hypothesis tests. At first the type I error proba- 
bility or significance level needs to be specified. This is the rate at 
which the test is allowed to falsely declare a data set as too dis- 
crepant to be compatible with the null hypothesis (the assumed dis- 
tribution function), even though it is in reality consistent. There is a 
value of the distribution of the test statistic, the critical value, which 
corresponds to this rate. If the value of the test statistic calculated 
from the data set then exceeds the critical value the null hypothesis 
is rejected for the data set. 

For some tests, the distribution of the test statistic can be cal- 
culated analytically for a fully specified null hypothesis, i.e. if no 
parameter is estimated. If parameters are estimated, the distribu- 
tion of the test statistic is not universal any more, but depends on 
the properties of the specific estimator. The distribution of the test 
statistic can then be obtained using a Monte-Carlo approach, of 
which follow the critical values. Typically the such derived critical 
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values are larger than for a fully specified hypothesis (cf. |Lilliefors| 
|1967| |1969| for the Kolmogorov-Smirnov test for the normal and 
exponential distribution). Therefore, if the critical values for the 
fully specified hypothesis are used when parameters are estimated, 
the results are conservative with an actually smaller type I error, but 
also less powerful. 

The significance level is not the only quantity characterising a 
statistical test. It can happen that a data set is not too discrepant to 
be rejected, but actually does not stem from the null hypothesis, i.e. 
a type II error occurs. The probability that a type II error does not 
occur is the (statistical) power of the test. If the power of the test 
is small then it is not very selective and the alternative hypothesis 
cannot be strongly excluded. For a given test the power can differ 
for various alternative hypotheses of the distribution. A demand for 
a general purpose test is to be powerful against a wide variety of 
alternative hypotheses. 

In order to be able to evaluate the "strength" of a statement 
concluded from a statistical test, it is therefore necessary to know 
the type I and type II error rates (the significance level and the 
power). However, not every test has the same power, therefore we 
conduct in what follows a power study which has the purpose of 
finding a powerful goodness-of-fit test to decide between infinite 
and truncated power laws. The astrophysical motivation for the 
choice of these hypotheses is the discussion in the literature about 
an upper mass limit for the distribution of stars in a star cluster (cf. 
|Weidner & Kroupa|2004||Oey & Clarke|2005||Koen|2006| > or about 
an upper limit for the star cluster luminosity function (cf. Gieles 
|et al.|2006| >. 



4.4 Description of the goodness-of-fit test statistics 

Goodness-of-fit tests can roughly be classified as tests based on 
the empirical DF (EDF), based on distance measures (e.g. the KS 
test) or the correlation coefficient, tests especially developed for 
a chosen null hypothesis (e.g. the Shapiro-Wilk test for exponen- 
tially), and tests to distinguish between two hypotheses (e.g. the 
Likelihood Ratio). Below we describe the test statistics used for a 
comparison. 

For the selection of the tests included in the comparison the 
properties of tests for exponentiality can be used, which can be 
found in the studies of Stephens i 1978 I, |D' Agostino & Stephens| 



(19861 and |Gan & Koehler| ( |1990hGan & Koehler| ( |1990| > which 
use EDF based tests included the alternative hypothesis of a trun- 
cated exponential, finding only very low powers. With the stabilis- 
ing transformation Kimber ( 1985 1 finds a larger power for the KS 
test, but did not include the truncated alternative. We include some 
EDF based tests in the original and stabilised version, as well as 
some tests based on tests for exponentiality, and two tests explic- 
itly for truncation. 



4.4. 1 EDF statistics based on distance measures 

The most prominent goodness-of-fit test is the Kolmogorov- 
Smirnov (KS) statistic (Stephens 1978; D'Agostino & Stephens 
[T986l|Gan & Koehler|1990| ~ 



D 



max 



i— 0.5 



1 

'2b' 



(20) 



which is the largest vertical distance between the data and the diag- 
onal in the PP plot (Rejection for D > D all ). The largest distance 



can also be measured in the stabilised PP plot, leading to the stabi- 
lized Kolmogorov-Smirnov statistic (Michael 1983, Kimber 1985), 



SD 



max 



i-0.5 



(21) 



|Kimber| {1985} found that for the exponential distribution these 
statistics are more powerful than the originals, but used a complete 
stabilising transformation (for both tails, So, eq. |18[l. Here only a 
right-tail-stabilising transformation (5, eq. \19\ is used, since it is 
more appropriate for the right-tailed power law and gives a better 
power. 

Another "measure of discrepancy" is the sum of the squared 
distances from the diagonal to the data point in a PP plot, the 
Cramer-von Mises statistic i Anderson & Darling|1952 ; Stephens 



[T978l|D' Agostino & Stephens|1986||Gan & Koehler|1990^ , 



C z 



E ho 

!=1 



2i-l 
In 



Yin 



(22) 



Like the KS statistic, this measure can be used in the stabilised PP 
plot, yielding the new stabilised Cramer-von Mises statistic, 

2 



sc 2 = £ (sp^-s 



2i- 1 
In 



(23) 



A modified form of the Cramer-von Mises statistic is the Anderson- 
Darling statistic, (Anders on & Darling||1952[ |Stephens||1978| 
|D'Agostino & Stephens|1986||Gan & Koehler|1990> 

A 2 = -L^ Z ^(log e (/ J (0 )-log f ,(l-P ( „ +1 _ /) ))-n, (24) 

!=1 " 

which gives more weight to the tails of the distribution. 



4.4.2 EDF statistics based on the correlation coefficient 
The correlation coefficient is a measure for linearity, given by 



(25) 



For perfect linearity R 2 has the value 1 or — 1 and for uncorrelated 
points R 2 = 0. If the points {^,y} are always positive as in our 
cases then R 2 lies in the interval [0,1]. The rejection criterion is 

R2 < R lnV 

The correlation coefficient can for example be used in the 
quantile-quantile plot, 



r 2 = R 2 [X^P {{) 



(26) 



but has, as the quantile-quantile plot, the disadvantage of needing 
the inverse distribution function. 

Another possibility is to use the correlation coefficient in the 
PP plot (PP correlation statistic ian & Koehler|1990) , 



r = R L I P {i) , 



i-0.5 



(27) 



or in the stabilised PP plot (stabilised PP correlation statistic, first 
proposed here), 



Sk z 



R 2 (SP (i) ,S 



i-0.5 



(28) 



The PP correlation statistic can be modified, as suggested by Gan 
& Koehler 1 1990 1, to force the points to go through {0.5, 0.5}. This 
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Table 1. Results of the power study. Tests are conducted under the null hypothesis of an infinite power law against the alternative hypothesis of a truncated 
power law with a type I error level of a/ = 0.05. The first column gives the size of the data set, n, and the second column the value of the exponent, a. The 
other columns give the power of the test statistic indicated in the top row in percent. The power is the fraction of the time in which the data drawn from the 
alternative hypothesis would be rejected (at the 5% level) as coming from the null hypothesis. The numbers on top of each group are the lower and upper limit 
of the parent distribution function. 



n 


J3 


D 




L 






,-2 


1-2 


1-2 






viz 


1 




A 


10. 150. 






























33 


1.7 


54.1 


72.0 


59.9 


69.4 


60.4 


51.7 


7.7 


59.9 


17.4 


65.9 


MA 


56.1 


64.6 


100.0 


50 


1.7 


70.9 


88.3 


75.0 


86.0 


78.7 


69.9 


5.1 


75.3 


16.3 


87.3 


68.7 


78.6 


86.3 


100.0 


99 


1.7 


93.8 


100.0 


96.8 


99.8 


98.6 


98.2 


3.7 


96.3 


40.3 


99.9 


96.7 


98.6 


99.9 


100.0 


33 


2.0 


23.3 


29.6 


24.4 


28.2 


26.0 


22.9 


4.5 


26.2 


7.0 


28.2 


22.2 


23.9 


29.8 


53.4 


50 


2.0 


25.0 


58.2 


32.8 


46.3 


36.2 


31.3 


3.7 


32.4 


7.8 


50.4 


38.8 


43.5 


61.0 


100.0 


99 


2.0 


39.1 


93.7 


50.7 


80.4 


63.7 


58.7 


3.1 


50.3 


11.5 


85.2 


68.9 


74.8 


94.9 


100.0 


33 


2.3 


7.8 


11.5 


8.9 


10.5 


9.2 


5.9 


5.4 


9.7 


5.5 


10.5 


10.7 


11.7 


11.3 


11.4 


50 


2.3 


8.5 


18.3 


8.8 


14.1 


9.3 


4.0 


2.9 


9.2 


3.7 


15.0 


14.2 


15.0 


18.9 


22.3 


99 


2.3 


11.9 


48.9 


11.5 


30.4 


15.6 


4.8 


4.2 


14.1 


4.9 


31.1 


29.3 


33.4 


53.5 


80.4 


10000. 1000000. 




























33 


1.7 


14.3 


18.5 


15.3 


19.1 


14.5 


18.8 


6.6 


15.3 


7.5 


19.3 


20.0 


17.9 


17.6 


19.5 


50 


1.7 


12.7 


33.3 


14.9 


25.6 


17.8 


28.3 


5.2 


18.4 


6.5 


26.5 


23.8 


25.6 


31.3 


50.9 


99 


1.7 


19.9 


77.4 


29.1 


57.3 


40.1 


60.1 


4.2 


30.9 


8.0 


61.8 


52.1 


58.0 


76.9 


100.0 


33 


2.0 


5.3 


5.7 


4.8 


5.6 


5.8 


5.9 


6.5 


4.8 


5.1 


5.9 


7.0 


6.8 


5.4 


6.2 


50 


2.0 


6.8 


8.4 


7.6 


7.8 


7.2 


5.5 


5.7 


8.3 


4.1 


7.3 


8.2 


8.8 


8.6 


8.7 


99 


2.0 


3.4 


14.5 


3.4 


6.2 


3.9 


5.3 


3.9 


3.7 


5.0 


9.8 


12.0 


11.7 


16.6 


12.6 


33 


2.3 


5.4 


4.9 


4.9 


4.4 


4.5 


3.3 


4.7 


4.1 


4.5 


4.2 


4.4 


4.6 


4.9 


4.9 


50 


2.3 


4.2 


6.0 


5.2 


5.9 


5.3 


2.2 


3.4 


4.4 


3.9 


5.5 


4.6 


4.9 


6.5 


5.5 


99 


2.3 


5.8 


7.0 


6.8 


6.9 


7.2 


1.1 


6.5 


6.0 


6.0 


6.9 


8.6 


8.6 


7.8 


6.8 



is done by repl acing in eq.|25|Z and Y by 0.5, denoting the modi 



fied version R 2 y 
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correlation statistic, 



(19901 found that the modified PP 
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is somewhat more powerful than k 2 . Again, the analogous pro- 
cedere is possible in the stabilised PP plot, giving the stabilised 
modified PP correlation statistic, 



Ski 



R o SP (i)> S 



i-0.5 



(30) 



4.4.3 Statistics based on tests for exponentiality 

Due to the connection of the power law and exponential distribu- 
tion, the various tests for exponentiality available in the literature 
are applicable for our purposes (cf. e.g. Be irlant et al.|2 006l. Since 
there is only a proportionality between p(x) and p e (log e x) most of 
the derived (exponential) null distributions for the following statis- 
tics are no longer valid for an infinite power law. 

The Shapiro-Wilk statistic {Shapiro & Wilk||1972| |Stephens] 
[T978l|D'Agostino&Stephens|1986| l, 



W 



-(X' 



-XL)Z(X!-X' 

i=l 



>\2 



(31) 



(X[ = log e Xi) is originally a two-sided statistic with minimum (n — 
I)- 2 and maximum 1 for the exponential case. For the use with a 
power law the rejection criterion for the alternative hypothesis of a 
truncated power law distribution is W > W cr i t in one-sided use. 
The Jackson statistic (|Jackson 1967| Stephens| |1978[ 



|D'Agostino & Stephens|1986||Beirlant et al.|2006 

Y" t- y" 

* ~ V" V" 



(32) 



with X" = log e (Xi/*Mi N ) and ?,-,„ = £^ =1 ^p^, is primarily the 
product of the ordered data and their expectation values XE(X^) = 
tj M , comparable to correlation type statistics. The division by 
Y!i=iXi removes the dependence on the scale parameter X. For a 
truncated power law alternative in one-sided use the rejection cri- 
terion is — T > — r cr ;t- 



Other tests for exponentiality are the the statistics of Brain & 
Shapir o] (|1983|> , the Moran statistic {Stephens|1978||D'Agostino & 
|Stephens|1986[ >, and the Greenwood statistic i Bartholomew! 1957} 
|stephens[l9*78l|D'Agostino & Stephens|1986) . We have also tested 
their powers when used for a power law, but they are not more 
powerful than the Shapiro-Wilk or Jackson statistic, and so we do 
not include details on them here. 



4.5 Tests for truncation 

The above described tests only allow one to distinguish whether 
the data are described by the null hypothesis or not. When the null 
hypothesis is rejected the test has to be made again, now with the 
alternative hypothesis as new null hypothesis. The Likelihood Ratio 
test combines this two-stage procedure into one test, by which can 
be decided whether of the two hypotheses is favourable. We use 
here the Likelihood ratio in the same way as the test statistics above. 
The test statistic is given by 



A = 



n"=iPpQ; 5 > £ MAx) 



(33) 



For the infinite case we used the ML estimate which does not in- 
clude a truncation to estimate the exponent and for the truncated 
case the modified ML for the exponent and upper limit. For numer- 
ical reasons the logarithm of eq.[33]is evaluated. 

An answer to the problem of estimating parameter and si- 
multaneously deciding between hypotheses can also be given in 
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the Bayesian framework of statistics which is, however, not in the 
scope of this article. 

A further specific test for truncation is the exceedance statistic, 



X 



max X, 



(34) 



(the largest data point), is only designed to test whether the distri- 
bution function is truncated or not. It cannot be used to detect a 
deviation from the power-law assumption. Furthermore it is one- 
sided with rejection criterion —X > — X^. 



4.6 Power comparisons 

The power of the various statistics was calculated at a signifi- 
cance/type I error level of (X/ = 0.05 with parameters for the power 
law as given in Table [T] The critical points for the null hypothe- 
sis of an infinite power law were calculated as follows. For each 
of the parameter combinations, but with xmax = °°, a Monte-Carlo 
sample containing 1000 data sets was generated. For each data set 
the parameters were estimated using the modified ML estimator 
(exponent eq. [T2] upper limit eq. |13} . Then the test statistics were 
calculated using the estimates when necessary. This gives the distri- 
bution of the respective test statistic, from which the critical value 
follows as the 95% quantile. 

For the power again a sample of 1000 data sets was generated, 
but now from a truncated power law. As before the parameters were 
estimated and the statistics calculated. The power is then the per- 
centage of data sets with a test statistic smaller than the critical 
value. 

The obtained powers are shown in Table [T] The exceedance 
statistic, X, is the most powerful test for truncation. However, it 
cannot be used for detecting deviations from the power law distri- 
bution. Thus it has to be used in conjunction with one or more of 
the other tests which include a test for the power law family as the 
parent distribution function. 

A general effect appearing for all statistics is that the power 
decreases with increasing slope and range of the limits. By such 
changes the truncated distribution becomes - informally speaking 
- more similar to the infinite distribution and thus harder to dis- 
criminate. Above a = 2 the performance of the tests drop signifi- 
cantly and therefore strong statements on truncation can barely be 
made. Unfortunately for the Salpeter value of the slope, a = 2.35, 
the studied tests are mostly not powerful enough to decide whether 
an upper truncation is present or not. However, in some not so ex- 
treme real cases such as the data set of massive stars in R136 a 
sufficient power can be achieved. Furthermore, even if a truncation 
cannot be detected then deviations from a power law might still be 
discoverable. 

Besides the general performance behaviour of the test statis- 
tics a further, rather surprising trend exists in the power. The most 
powerful tests are not necessarily the tests derived especially for the 
power law distribution from tests for exponentiality. The stabilising 
transformation (eq. |19} strongly enhances the power of general- 
purpose ECDF or correlation statistics so that they outperform 
the specialised tests. The Kolmogorov-Smirnov statistic which is 
known to be not very powerful (cf. Gan & Koehler 19901 becomes, 
after stabilisation, more powerful than all other tests except for the 
exceedance test. In their not stabilised forms the general-purpose 
tests are, as expected, less powerful than the specialised tests. This 
enhanced power is a useful property since general-purpose statis- 
tics can easily be modified to tests for a different null hypothesis, 
e.g. a two-part power law. 



Table 2. Estimates for the 29 most massive stars in R136. The standard 
deviations of the estimators were calculated using a Monte-Carlo sample of 
size 10 000. For the binning methods 5 bins are used. The bias was calcu- 
lated using the results of the modified ML method as input values. 

Data using Chlebowski & Garmany 1 1991 1 



Estimate for 


Slope 5 


Bias slope 


Mmax 


Bias Mmax 










[M e ] 




Const. Bins, LR 


3.38±0.72 


0.43 






Var. Bins, y} 


2.42±0.75 


-0.01 


134±12 


-11 


CCDF 




2.02±0.88 


-0.11 


140±9 


-1.9 


Beg 




2.17±0.77 


<0.01 


142±8 


-0.4 


Beg, recursive 


2.17±0.77 


0.01 


142±8 


-0.4 


ML°° 




3.51±0.35 


1.34 






ML 




2.11 ±0.73 


-0.06 


136±7 


-8 


Mod. ML 




2.20±0.78 


0.02 


143±9 


< 0.1 


Results of 


Koen 


{2006} 








CCDF 








143.9 




ML 




2.11 




136 





In summary, the best test for truncation is the exceedance test, 
X. To confirm the hypothesis of a power law and for better signif- 
icance this test should be followed by some of the most powerful 
remaining tests. These are, loosely ordered in descending power, 
the stabilised Kolmogorov-Smirnov test SD, the stabilised PP cor- 
relation test S&q, the stabilised Cramer- von Mises test SC 2 and the 
Jackson statistic T. 

When a truncation is detected, then the hypothesis of a trun- 
cated power law has to be confirmed by again applying the respec- 
tive statistics with this distribution (the truncated power law) as the 
null hypothesis. 



5 EXAMPLES 

5.1 The massive stars in R136 

As a first exemplary application of the presented statistical tech- 
niques, in particular of the estimators, we chose the data set of 
massive stars in R136 published by Ma ssey & Hunter] ( |1998| >. 
They gave for the 29 most massive stars the masses based on 
two different stellar models ( jChlebowski & Garmany||1991| with 
masses ranging from 56 Mq to 136 M Q , and |Vacca et al. [l996[ 
75 — 155 Ma). The results of the estimators are shown in Table 
where a Monte-Carlo sample of size 10000 was used to calculate 
the standard deviations. 

Beg's estimator and the modified ML method agree well (S = 
2.2), the ML estimate is slightly smaller (5 = 2.1). In reasonable 
agreement with this value are also the results of variable-size bin- 
ning and fitting the complementary cumulative DF plot. For com- 
parison the results of Koen ( 2006 ) are also given in Table [2] The 
ML estimates are equal, only the CCDF result differs, likely due 
to a different definition of the empirical DF (Koen uses i/(n+ 1) 
whereas here (('— 0.5)/« is used). The ML method without includ- 
ing an upper limit gives a much larger exponent (S = 3.5) which 
shows the effect of a model mismatch. A comparison of this value 
only with a constant-size histogram, where a linear regression gives 
S = 3.4, would not give any indication of the mismatch. 

The upper limit is determined as ~ 140 M s by Beg's estima- 
tor, the ML and the CCDF method. The results from variable-size 
binning are not consistent with the data set, because this upper limit 
is smaller than the largest data point. 
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Figure 6. Truncated SPP plot of the massive stars in R136 with masses 
according to the model of Chlebowski & Garmany ( 1991 1 and parameters 
estimated using the modified ML method. Also shown is the curve for a in- 
finite power law (dotted). The parallels to the diagonal limit the acceptance 
region of the stabilised KS test, null hypothesis of a truncated power law, 
significance level 5%. 



to 
o 

3 s 
to 



logio M e 
2.5 3.0 3.5 4.0 4.5 




S(PUM (i) j) or S(P(M)) 



Figure 7. Infinite SPP plot of the LMC star clust ers (age < 10 7 5 yr) with 
the lower mass limit of |de Grijs & Anders|{2006| , 10 2 2 M s . (Dotted line: 
infinite hypothesis, 2mml = 1.47; solid curve: truncated hypothesis, pa- 
rameters as estimated (a = 1.47, Mmax.mml — 64200 Mp ; ); dashed curve: 
truncated hypothesis, a = 2, Mmax.mml — 64200 M,: : ; dash-dotted lines: 
limits of the acceptance region of the stabilised Kolmogorov-Smirnov test, 
significance level 5%). 



For the goodness-of-fit analysis an SPP plot with a truncated 
power law as null hypothesis is shown in Fig. [6] The curve for 
the infinite power law is clearly not fitting the data. The stabilised 
Kolmogorov-Smirnov, Cramer-von Mises and SPP correlation co- 
efficient test all give a strong disagreement of the data with an infi- 
nite power law and no disagreement with a truncated power law. 

The mod ified ML estimates from the data set using the models 
of |Vacca et al.|l996| are 5 = 2.87 ±0.98 and M M AX = 163±9M Q . 
The goodness of fit tests indicate a truncated power law with high 
significance too. 

5.2 The young star clusters in the Large Magellanic Cloud 

The second example for the methods presented above, with an em- 
phasis on the advantage of the SPP plot, is the analysis of the mass 
distribution of young star clusters in the Large Magellanic Cloud. 
We use a part of the data set given by de Grijs & Anders (20061, 
the star clusters with ages younger than lO'' 5 yr and masses larger 
than 10 2 ' 2 M Q . 

Based on an inspection of the shape of a histogram of the data 
|de Grijs &"A nders (2006) concluded that there is a significant flat- 
tening of the mass function for M < 10 3 Mq (see their fig. 8). In- 
deed, also an SPP plot with an infinite null hypothesis, Fig|7] shows 
that the empirical curve of the data is strongly bent in the lower 
mass range (M < 10 2 ' 5 Mq, estimated exponent ~ 1.5). In the up- 
per mass range the infinite SPP plot reveals that the data are better 
described by a truncated power law (solid line). This indicates that 
above ~ 10 2,5 M K , the data presumably will be consistent with a 
truncated power law. 

With an SPP plot using only the star clusters more massive 
than ~ 10 2,5 Mq this hypothesis is confirmed (Fig. [5|, the sta- 
bilised Kolmogorov-Smirnov acceptance region is not exceeded. 
Thus, a change in the slope or shape of the mass function in the 
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Figure 8. Truncated SPP plot of the LMC star clusters (age < 10 7 5 
yr) starting at 10 25 M©. (Dotted line: infinite hypothesis, 0(mml = 
1.62; solid line: truncated hypothesis, parameters as estimated (Smml = 
1.62, Mmax.mml = 68000 Mg ); dashed: truncated hypothesis, a = 2, 
Mmax.mml = 68000 M ; dash-dotted: limits of the acceptance region of 
the stabilised Kolmogorov-Smirnov test, significance level 5%). 
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Figure 9. Influence of the completeness (dotted line) on the observable 
mass function (solid line), based on an assumed power law (dashed line) as 
the underlying distribution function. The parameters were chosen to match 
the situation for the Large Magellanic Cloud, see text. The vertical lines 
correspond to a: 10 1 58 M a , b: 10 1 98 M Q , c: I ()" ' M a , and d: 10 2 5 M . . 
The mass function is scaled arbitrarily for better visibility. 



mass range 

10 2.5 _ ,q3 Mq> as stated by |de Grijs & And ers (20061, 
cannot be deduced using our techniques. Only the mass range 
10 2 - 2 — 10 2 ' 5 Mq seems to deviate from the power law. The slope 
which is derived from the data with masses larger than 10 2 ' 5 M 
is a = 1 .6 ± 0.06 and an upper limit 68 ± 1 1 x 1O 3 M is obtained 
by the modified Maximum Likelihood method. This exponent is 
smaller than the value determined by de Grijs & Anders (20061, 
a — 1.8 ±0.1, who used constant-size binning and star clusters 
more massive than 10 3 Mg (for this mass range the modified max- 
imum Likelihood estimate is a = 1 .63 ±0.1). 
The feature in the mass range 

10 2.2 _ 10 2.5 Me could be 

caused by an actual change of the mass function. However, since 
it is at the lower mass end, it could also be caused by an incomplete 
data set. The completeness limit adopted by |de Grijs & Anders] 
(2006 ) was derived by jHunter et al. l ^2003} from the b ehaviour of 
the luminosity function (see fig. 4 of Hunter et al. 2003). They used 
as the brightness limit the brightness where the luminosity function 
reaches at the faint side half of its peak value, obtaining ^My = 
-3.5 mag or 10 L58 M (using M = l0 6 +°- 4 (- 14 - 55 --^) m q< 
|Hunter et al.|2003| eq. 1). The mass related to the brightness limit 
is valid for clusters of an age of 10 Myr. In a similar way |Parmen-| 
|tier & de Gri js ( 2008 1 derive - starting from the mass distribution 
of a chosen age interval older than 10 Myr - from the mass which 
separates the lower 25% from the upper 75% a completeness limit 
of = —4.7 mag. If we use this value also for younger clusters, 
then a completeness mass of 10 2 06 M would result (However, for 
unknown reasons an application of their method to clusters younger 
than 10 Myr leads to a different completeness mass of 10 2 - 35 M©, 
Parmentier, priv. comm.). A completeness limit derived in such a 
way coincides approximately with the peak of the observed mass 
function. But the transition from no detection to complete detec- 
tion is smooth and has a certain broadness in which only a fraction 
of all sources is detected which can affect a wide mass range, il- 
lustrated in Fig. [5] The observable mass function (solid line) is the 
product of the actual mass function (dashed, exponent a = 1.6) and 
the completeness function (dotted). As the functional form for the 



c(M) 



M 



(35) 



The parameters of the completeness function were chosen such that 
the half peak point of the observable mass function is at 10 158 M@ 
(or ^-y = — 3.5mag, as Hunter et al. 2003 . point a in Fig_j9j and the 
peak mass is ~ 10 2 Mq (^V ~ —4.5, point b in Fig.lpJTlt is just 
a coincidence that for the used parameters (log IO Mo = 1.98 and 
(j> = 3.12) the 50% completeness mass of the completeness func- 
tion coincides with the peak mass of the observable mass function. 
With these empirically determined parameters the observable mass 
function is shallower than the actual power law in the mass range 
below fs 10 2 ' 5 Mq. This strongly supports the argument that the 
deviation of the data from the power law in Fig.|7]is caused by in- 
completeness. The distribution of star clusters with ages < 10 yr 
are well consistent with a single power law with 2 = 1.6, starting 



6 SUMMARY AND CONCLUSIONS 

In this work we compared methods to estimate the exponent 
and upper limit of a truncated power law distribution. We re- 
viewed graphical methods to represent the data. Finally we studied 
goodness-of-fit tests, specifically to test for truncation. 
Our results are: 

(i) A generally working estimator for the exponent and upper 
limit is our modified maximum likelihood method. It performs well 
with respect to bias and standard deviation. 

(ii) A maximum likelihood estimate of the exponent without 
considering a truncation can lead to biased results if the data stem 
from a truncated power law. 

(iii) The estimator of Beg ( 1983 1 is also performing well but is 
numerically not stabl e. Variable-size binning as introduced by |Mafz| 



Apellaniz & Ubeda 1 2005 i performs well for the exponent. The 
estimate for the upper limit based on the normalisation constant is 
biased. 

(iv) The stabilising transformation introduced by |Michael| 
( 1983) enhances plots and goodness-of-fit tests. For one-sided dis- 
tributions only a half transformation should be made to achieve op- 
timal results. 

(v) The stabilised PP plot is a particular useful display of the 
data. 

(vi) The stabilised Kolmogorov-Smirnov statistic (SD), the sta- 
bilised PP correlation test (Sk 2 ), the stabilised Cramer-von Mises 
statistic (SC 2 ), the Jackson statistic (T) and the QQ correlation (r 2 ) 
test are powerful goodness-of-fit tests for the truncated power law. 

(vii) The exceedance statistic (X) is the most powerful test for 
truncation. Since it does not test for power-law behaviour it has to 
be used in combination with a powerful goodness-of-fit test for the 
truncated power law, as the ones mentioned in the previous point. 

(viii) The massive stars in R136 are well described by a trun- 
cated power law with S = 2.20 ± 0.78 and M M AX = 143 ± 9 M , 
using the Chlebowski & Garmany (1991) stellar models for mass 
determination, or a = 2.87±0.98 andM M AX = 163±9 M , using 
the |Vacca et al.| (T996l stellar models. 

(ix) The young star clusters in the Large Magellanic Cloud (ages 
younger than 10 7,5 yr) with masses larger than 10 2,5 M are well 
described by a truncated power law with S = 1 .62 ± 0.06 and 
M MAX = 68.8 ±11.6 xlO 3 M . 



1 2 Th. Maschberger and P. Kroupa 



(x) A change in shape of the star cluster mass function in the 
Large Magellanic cloud in the low mass range M < 10 3 M@, as 
reported by |de Grijs & An ders ( 2006), cannot be verified. For 
M > 10 2,5 Mq the observed distribution follows a truncated power 
law, a flattening below 10 2,5 Mq is most likely caused by an under- 
estimated completeness limit. 



Vacca W. D., Garmany C. D., Shull J. M., 1996, ApJ, 460, 914 
Weidner C, Kroupa P., 2004, MNRAS, 348, 187 
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7 ACKNOWLEDGEMENTS 

We thank Cathie Clarke and Douglas Heggie for critical reading of 
the manuscript and valuable comments. TM acknowledges finan- 
cial support by the AlfA. 



The estimator for the exponent (Beg 1983) is given in its original 
form as 



(«-3)!i;.: () (-iyr] 2 K;i 4 



(»-4)!Ef =0 (-iy( n 7 2 )^r 



(Al) 
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where = S - 1 and Kj = T - nY - j(Z - Y) with Y = log e X (1) , 
Z = log e Xj„) and T = £" =1 log e Xj. The terminating index of the 
sum, j*, is determined by the condition T — nY — j(Z — Y) > 0, as 
shown by Beg 1 1983 i. The estimate for the exponent is then S = 
+ 1. 

The direct evaluation of eq. |A1| involves the calculation of 
("J 2 ) , which is only practicable for less than about 170 data points 
in double precision arithmetic. This problem can be handled with a 
recursive implementation of the estimator, feasible for any number 
of data, as follows. 



To abbreviate we introduce L 



V 1 J T-nY , 



K' 



n—A 



which leads to 
(A2) 



With changing the limits of the sum and omitting (n — 3)! the nu- 
merator of eq. |Al| reads 



(r-«yr 4 £(-i)^J]i^ 
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(A3) 



Omitting the prefactor (T — nY) n 4 , the expanded sum reads 



(n-2)! 
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Starting with the second term this can be written as 



(A4) 
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(A5) 



The superscript (n — 4) should only indicate the exponent and is not 

used as an exponent in Sj 4 ' . From this the recursion can easily 
be seen: 



c("-4) 
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where / descends from j* to 2 and 5 ; > = 0. The last step is 
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(A7) 



The recursion for the denominator in eq. |Al| is as for the numerator, 
but replacing the exponent n — 4 by n — 3 in equations |A6| and |A7| 
The estimator of is then (remembering all omitted factors) 



n-3 S("- 4 > 
T-nY S("" 3 ) ' 



(A8) 



Estimators and goodness-of-fit tests for power 



The estimators for the upper limit in the form of Beg ( 1983 1 is 

1 i ij: (-ivcy)*r 2X 
, ( n - i )E^ (-iy("7 2 )^? 



The recursion formula for the sum in the numerator follows by 
analogous steps as before with 

= n —l(L^-s\ n - 2) ), (A10) 

the last step 

S '(n-2) = 1 _ 5 >-2). (AU) 

and 

re(n-l) S("- 3 ) 
with S'"~ 3 ' from the estimator for the exponent above. 



-max = X {n) (l+ n(n _ 1} 5(n _ 3) ), (A12) 



This paper has been typeset from a TjjX/ IATgX file prepared by the 
author. 



